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Fair  Polyline  Networks  for  Constrained  Smoothing 
of  Digital  Terrain  Elevation  Data 

Michael  Hofer,  Guillermo  Sapiro,  Senior  Member,  IEEE,  and  Johannes  Wallner 


Abstract — We  present  a  framework  for  smoothing  grid-like 
digital  terrain  elevation  data,  which  achieves  fair  shape  by 
means  of  minimizing  an  energy  functional.  The  minimization 
is  performed  under  the  side-condition  of  hard  constraints  which 
come  from  available  horizontal  and  vertical  accuracy  bounds 
in  the  elevation  specification.  We  introduce  the  framework  and 
demonstrate  the  suitability  of  this  method  for  the  tasks  of 
accuracy-constrained  smoothing,  feature-preserving  smoothing, 
and  filling  of  data  voids. 

Index  Terms — Digital  terrain  elevation  data  (DTED),  surface 
smoothing,  fair  polyline  networks,  guaranteed  error  bounds, 
topography  preserving. 


1.  Introduction 

SMOOTHING  digital  terrain  elevation  data  (DTED)  with 
guaranteed  error  bounds  and  feature  preservation  is  of 
great  importance  in  practical  geoscience  tasks,  both  for  visu¬ 
alization  purposes  and  for  post-processing  (e.g.,  compression). 

A  typical  data  set  we  are  working  with  is  depicted  in  Fig.  1. 
We  see  that  the  iso-height  contours  are  often  jagged,  and  false 
contours  appear  (small  “islands”).  This  is  due  in  part  to  the 
high  level  of  detail  present  in  the  data  and  in  part  to  data 
acquisition  and  generation  limitations.  The  noisy  aspect  of  the 
DTED  can  also  be  seen  by  means  of  a  shaded  relief  (Fig.  2). 
After  constrained-smoothing  with  the  technique  here  proposed, 
iso-height  contour  lines  look  like  those  in  Fig.  3. 

The  data  we  use  have  been  provided  by  the  U.S.  National 
Geospatial-Intelligence  Agency  (NGA).  The  data  have  been 
obtained  by  the  Shuttle  Radar  Topography  Mission  (SRTM) 
in  February  2000,  one  of  the  most  significant  space  surveys 
of  earth  ever  undertaken,  see  e.g.,  http ://www2.jpLnasa. gov/ 
srtm/  Sind  [1].  SRTM  is  a  joint  project  between  NASA,  NGA, 
the  German  Aerospace  Center  (DFR),  and  the  Italian  Space 
Agency  (ASI).  The  data  complies  to  the  DTED-2  specification, 
i.e.,  is  given  as  a  uniform  gridded  matrix  of  terrain  elevation 
values  with  a  spacing  of  one  arc  second  (approximately  30 
meters).  The  specification  also  includes  information  about 
absolute  errors  (or  accuracy),  both  in  horizontal  as  well  as  in 
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Fig.  1.  Iso-height  contour  plot  of  original  elevation  data. 


Fig.  2.  Shaded  relief  of  original  data. 


Fig.  3.  Iso-height  contour  plot  of  smoothed  elevation  data. 

vertical  directions.^  This  means  that  the  data  includes  specific 
information  about  the  possible  errors  in  the  posts  positions 
(circular  horizontal  error)  and  the  reported  height  (vertical 
error).  This  additional  accuracy  information  is  exploited  by 
our  proposed  smoothing  framework. 

^For  details  on  the  specification,  see  for  example  http://www.fas.org/ 
irp/program/core/dted.  htm.  See  also  http: //mountains,  ece.  umn.  edu/r^ guide/ 
dtedspecification.pdf  for  a  copy  of  the  official  DTED  specification. 
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It  is  of  great  importance  that  smoothing  digital  elevation 
data  respects  the  given  error  bounds.  In  other  words,  the 
smoothed  data  must  be  in  between  the  absolute  accuracy 
bounds  provided  by  the  specification.  This  is  needed  in  order 
to  guarantee  that  the  processed  data  is  in  the  allowed  (by  the 
data  specification)  height  and  position  range.  In  this  paper 
we  present  an  algorithm  which  does  just  that,  and  is  also 
capable  of  filling  holes  as  well  as  preserving  features  in 
a  consistent  manner.  The  hard  error  (or  accuracy)  bounds 
are  represented  as  tolerance  cylinders,  where  the  surface  is 
permitted  to  move.  The  smoothing  is  accomplished  by  means 
of  an  energy  minimization  technique  (fair  polyline  network), 
constrained  to  keep  the  surface  inside  the  tolerance  bounds. 
Relaxing  or  limiting  the  size  of  the  cylinders,  void  filling  and 
feature  preservation  are  naturally  obtained  within  the  same 
framework.  The  proposed  technique  can  be  applied  both  to 
smooth  the  whole  data  (as  a  function  on  the  plane)  or  singular 
iso-height  lines. 

The  remainder  of  the  paper  is  organized  as  follows.  In 
Section  II  we  briefly  describe  related  prior  art.  Our  contri¬ 
butions  are  put  forward  in  Section  III.  In  Section  IV,  we  first 
define  fair  polyline  networks  and  then  discuss  the  error  bounds 
and  the  smoothing  procedure.  For  validation  of  the  algorithm, 
we  introduce  a  deviation  measure  based  on  a  gauge  body. 
Furthermore,  we  show  that  besides  smoothing,  our  algorithm 
is  also  capable  of  preserving  features  and  filling  voids.  In 
Section  V  we  present  experimental  results,  and  in  Section  VI 
we  conclude  the  paper. 

II.  Previous  work 

The  present  paper  is  motivated  by  our  work  on  variational 
interpolation  of  subsets  [2],  on  energy-minimizing  splines  in 
manifolds  [3],  and  on  fair  curve  and  polyline  networks  in  non¬ 
linear  geometries  [4],  [5].  In  the  following  we  review  literature 
related  to  surface  smoothing,  guaranteed  error  bounds,  and 
hole  filling,  with  a  focus  mainly  on  the  geometry  processing 
rather  than  on  the  geosciences  community. 

A.  Surface  smoothing 

For  data  of  a  surface  like  terrain,  the  filtering  of  noise 
and  smoothing  of  the  geometry  has  been  of  great  interest. 
Constraints  imposed  on  the  smoothing  process  include  the 
preservation  of  linear  and  non-linear  surface  features  such 
as  sharp  edges,  corners,  or  non-planar  curves.  This  can  be 
achieved  for  example  using  the  so-called  anisotropic  smooth¬ 
ing  methods  (cf.  [6]  and  [7])  in  this  area.  There  are  only  few 
publications  where  hard  error  bounds  are  incorporated  into 
these  geometric  regularization  methods.  Recently,  geometric 
active  contours  have  been  used  as  the  base  of  a  constrained 
regularization  framework  for  digital  elevation  data  [8]. 

B.  Error  bounds  and  the  tolerance  aspect 

While  there  are  many  contributions  concerning  error  propa¬ 
gation  in  digital  elevation  models  (see  for  example  [9]  and  the 
references  therein),  it  is  apparently  difficult  to  locate  previous 
work  on  smoothing  of  terrain  data  with  guaranteed  error 


bounds  (being  the  only  paper  to  the  best  of  our  knowledge, 
our  work  in  [8]).  Thus,  it  is  perhaps  no  coincidence  that  the 
present  contribution  to  this  topic  has  its  origin  in  geometry 
processing  and  geometric  modeling.  We  would  also  like  to 
mention  that  the  topic  is  related  to  tolerance  analysis:  It  is 
common  to  locate  imprecisely  defined  entities  in  computations 
by  tolerance  zones.  The  most  prominent  example  of  this  is 
to  compute  with  intervals  instead  of  with  real  numbers  [10], 
[11].  Geometric  operations  with  tolerance  zones  have  been 
studied  later  [12].  In  that  paper,  energy-minimizing  curves 
which  interpolate  imprecisely  defined  points  (i.e.,  tolerance 
zones)  are  considered,  which  is  also  the  topic  of  [2],  and  which 
is  closely  related  to  the  problem  solved  in  the  present  paper. 

In  the  approximation  theory  literature,  the  problem  of 
regularization  with  constraints  has  been  addressed  by  Kimel- 
dorf  and  Wahba  [13].  They  showed  how  to  compute  one 
dimensional  splines  with  hard  vertical-error  constraints.  While 
this  elegant  approach  can  easily  be  extended  to  higher  di¬ 
mensions,  it  doesn’t  include  the  horizontal  freedom  given 
by  the  horizontal  absolute  error.  It  is  also  not  developed  for 
the  additional  geometric  constraints  that  are  natural  to  add 
in  our  framework.  The  theory  of  total  least  squares  [14] 
also  addresses  the  “freedom  of  motion”  of  the  given  data, 
both  in  the  vertical  and  horizontal  position.  In  it’s  original 
form,  although  computationally  very  efficient,  the  framework 
does  not  provide  hard  constraints  (that  is,  the  error  is  not 
guaranteed  to  be  below  the  allowed  bounds),  neither  does  it 
include  any  kind  of  explicit  regularization  or  geometric  terms. 
In  order  to  add  these  important  constraints,  the  problem  has 
to  be  transformed  into  a  variational  formulation,  much  of  the 
flavor  here  introduced.  The  recent  work  in  [15]  presents  the 
problem  of  level-set  estimation  as  a  tree  optimization  one,  with 
guaranteed  optimality  but  no  guaranteed  error  bounds. 

C.  Void  filling 

It  is  possible  that  in  some  areas,  no  elevation  values  are 
available  -  we  refer  to  these  void  areas  as  ‘holes.’  There  are 
various  reasons  for  holes  in  the  data.  They  can  be  caused,  e.g., 
by  occlusion  or  poor  signal  returns. 

The  literature  on  hole  filling  algorithms  for  geometric  mod¬ 
els  roughly  fits  in  three  categories:  surface  based,  volumetric, 
and  example  based  methods.  The  article  [16]  surveys  the 
literature  on  the  first  two  approaches  up  to  the  year  2002.  Ideas 
from  image  inpainting  are  used  in  the  volumetric  approach  by 

[17] .  Filling  holes  in  point  set  surfaces  is  also  discussed  in 

[18] .  Recently,  example  based  methods  have  emerged  [19]- 
[21].  It  is  interesting  to  note  that  the  example  based  approach 
is  most  closely  related  to  a  common  practice  of  filling  voids 
in  elevation  data:  One  uses  data  from  alternate  sources  and 
blends  them  with  the  data  surrounding  the  void  (the  ‘fill  & 
feather  method,’  cf.  [22],  [23]). 

III.  Contributions  of  the  present  paper 

Elevation  data  of  the  DTED-2  specification  are  given  as  a 
height  field  over  a  uniform  grid  in  the  x^-plane,  such  that 
every  data  point  has  horizontal  and  vertical  error  (accuracy) 
bounds.  Eig.  4  shows  data  points  together  with  their  error 
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Fig.  4.  An  initial  surface  interpolating  the  given  data  points  The  surface 
can  be  smoothed  such  that  it  stays  within  all  tolerance  cylinders  Zi^  of 
radii  and  height  2hij.  The  radii  and  height  are  specified  with  the  DTED 
information.  A  different  size  of  the  error  cylinders  can  be  used  to  preserve 
surface  features  (e.g.,  ridges)  during  the  smoothing. 

bounds  in  the  form  of  cylinders  of  revolution  Zij ,  and  an  initial 
surface  that  interpolates  the  given  points  The  objective  is 
to  compute  a  smooth  surface  that  stays  within  the  given  error 
margins  dictated  by  the  cylinders. 

The  framework  we  propose  views  the  height  field  over 
the  x^-plane  as  a  network  consisting  of  x-parallel  and  y- 
parallel  polylines,  as  shown  by  Fig.  5.  We  perform  smoothing 
by  minimizing  a  discrete  bending  energy  of  these  polylines, 
always  respecting  the  error  bounds.  Our  framework  is  very 
flexible  and  allows  e.g.,  to  use  individual  error  bounds  for 
every  data  point,  fill  voids  in  the  data,  or  smooth  only  a  subset 
of  the  data  for  feature  preservation.  We  later  illustrate  this  at 
hand  of  several  examples. 

IV.  Fair  polyline  networks  eor  terrain  smoothing 

WITH  GUARANTEED  ERROR  BOUNDS 

We  first  define  fair  polyline  networks  and  then  we  show 
how  to  extend  them  for  smoothing  of  digital  terrain  elevation 
data  with  guaranteed  error  bounds. 

A.  Fair  polyline  networks 

Our  smoothing  procedure  is  based  on  minimizing  the  energy 
of  a  polyline  network.  We  first  define  the  energy  of  a  single 
polyline,  and  then  the  energy  of  the  polyline  network  as  the 
sum  of  energies  of  all  the  polylines  that  contribute  to  the 
polyline  network.  For  a  smooth  curve  c(t),  defined  in  some 
parameter  interval  [a,  6],  the  linearized  bending  energy  (cubic 
spline  energy)  is  defined  by 

E{c)=!l\\c{t)rdt.  (1) 

Its  minimizers  under  interpolation  conditions  are  the  cubic 
B-spline  curves  [24].  A  polyline  p  =  (qi,q2,  •  •  •  ,qL),  as  a 
discrete  curve,  possesses  a  discrete  linearized  bending  energy: 

-E  =  IIA^qjlP,  A^qi  =qj_i -2qj  +  qi+i.  (2) 

Curve  networks  and  polyline  networks,  which  are  the  topic 
of  [4],  [5],  have  energies  which  are  defined  as  the  sum  of 
the  energies  of  the  curves  or  polylines  which  they  are  made 


Eig.  5.  A  polyline  network  with  vertices  pij  and  two  families  of  polylines 
(solid  and  dotted)  representing  a  height  field. 

of.  The  given  elevation  data  constitute  a  rectangular  array  of 
points:  Zij),  (i  =  j  = 

We  define  the  energy  of  the  data  collection  to  be  the  sum  of 
energies  of  the  N  different  polylines  defined  by  j  =const.  and 
the  M  different  polylines  defined  by  i  =const.: 

'^zZ  -^zZ  -o  IIPlj-1  “  +  Plj+iII  • 

1 — r  J — z 

A  fair  polyline  network  is  one  which  has  minimal  energy 
among  all  polyline  networks  which  fulfill  a  fixed  set  of 
constraints. 

B.  Error  bounds 

The  data  sets  we  use  consist  of  points  p^y  which  are 
equally  spaced  in  x  and  y  direction.  We  might,  without  loss  of 
generality,  assume  that  Xij  =  g  -  i  and  yij  =  g  '  j,  where  g  is 
the  grid  element  size  (approximately  30  meters  for  DTED- 
2).  The  data  points  contain  errors,  both  in  horizontal  and 
vertical  direction.  If  maximum  horizontal  error  of  a  data  point 
Y>ij  is  bounded  by  and  the  vertical  error  by  hij,  then 
the  true  location  of  that  data  point  is  within  a  cylinder  Zij 
of  diameter  2rij  and  height  which  is  centered  in  the 

given  point  p^y  =  {xij ,  yij ,  Zij ) }  Each  point  p^ y  is  equipped 
with  its  own  tolerance  cylinder  Zij  (Fig.  6,  left).  We  also 
refer  to  these  tolerance  cylinders  as  error  cylinders.  They  are 
hard  constraints,  which  means  that  the  terrain  surface  we  are 
seeking  has  to  pass  through  all  Z^y’s. 

Fig.  6,  top  right,  shows  the  initial  state  of  the  network,  only 
one  polyline  is  shown  as  a  representative  of  the  surface  defined 
by  the  data.  Our  goal  ultimately  is  to  move  the  points  p^y 
such  that  the  energy  (3)  becomes  smaller,  but  the  surface  still 
passes  through  the  cylinders  Zij.  We  can  achieve  this,  e.g., 
by  requiring  that  the  points  p^y  themselves  do  not  leave  Z^y, 
as  illustrated  by  Fig.  6,  at  right,  center.  This  condition,  which 
is  referred  to  as  Option  7,  however,  is  stricter  than  actually 
necessary.  It  is  certainly  sufficient  that  for  each  vertex  p^y, 
one  of  the  two  polylines  meeting  there  meets  Zij  (Fig.  6, 
bottom  right).  This  slightly  more  relaxed  condition  is  referred 
to  as  Option  2. 

^Other  error/accuracy  models  will  just  define  different  tolerance  shapes, 
without  altering  the  framework  here  presented.  Cylinders  are  the  appropriate 
shape  to  represent  the  absolute  errors  specified  with  DTED. 
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Fig.  6.  Smoothing  height  fields  via  fair  polyline  networks.  (Left)  Polyline 
network  with  the  tolerance  cylinders  Zij  associated  with  the  vertices  pij. 
(Right)  At  top,  initial  state  before  optimization.  At  center,  Option  I  of 
optimization  (all  vertices  stay  inside  the  tolerance  cylinders).  At  bottom, 
Option  2:  The  vertices  are  allowed  to  leave  the  cylinders,  provided  the  surface 
still  passes  through  them. 

C.  The  smoothing  procedure 

Minimizing  the  quadratic  function  (3)  of  the  variables 
Xij^Pij^Zij,  subject  to  the  constraints  mentioned  above,  is  a 
quadratic  programming  problem  with  non-convex  side  condi¬ 
tions.  It  contains  too  many  variables  for  just  submitting  it  to 
a  generic  optimization  procedure  as,  e.g.,  provided  by  math¬ 
ematical  software.  The  following  properties  however  allow  a 
more  direct  approach: 

(a)  If  computed  with  x  and  y  coordinates  of  the  initial 
data  alone,  the  energy  would  be  zero,  owing  to  the  regular 
alignment  of  initial  data.  We  therefore  don’t  expect  the  points 
Pij  to  move  very  much  in  x  or  ^  direction  during  optimization. 
Indeed  this  is  confirmed  by  numerical  experiments. 

(b)  At  each  location,  the  smooth,  energy-minimizing  terrain, 
is  defined  by  the  input  data  which  are  nearby.  Thus  there  is  no 
need  for  globally  minimizing  the  energy  over  large  data  sets. 

(c)  The  condition  that  pij  remains  inside  Zij  (Option  1) 
leads  to  a  convex  optimization  problem  which  has  a  unique 
solution.  As  illustrated  by  Fig.  6  (right),  the  difference  between 
Option  1  and  Option  2  is  not  very  big.  This  means  that  the 
optimization  problem  we  have  to  solve  when  using  Option  2  is 
‘convex  enough’  so  that  we  don’t  expect  a  direct  minimization 
procedure  getting  stuck  in  a  local  minimum. 

In  view  of  property  (a),  we  only  use  the  z  coordinates  of 
the  data  points  as  variables  for  minimization: 

J  —  i  %  —  z 

^  ~  •  (4) 

Because  of  (c),  we  employ  a  gradient  descent  method,  with 
the  original  elevation  data  as  initial  condition.  It  is  elementary 
that  the  gradient  of  the  energy  is  given  by 

(S/E)ij  =  {Zi-2J  +  Zij-2,j  +  ^i,j-2  +  ^i,i+2) 

provided  j  >  2  and  i  <  M  —  I,  j  <  N  —  1.  For  vertices 
near  the  boundary  there  are  different  formulae. 

Optimization  is  basically  implemented  as  follows:  First,  we 
find  a  direction  of  descent,  e.g.,  by  letting  pij  := 


We  consider  the  1 -parameter  variation  Ef  of  the  energy  (4) 
defined  by  the  z  coordinates  Zij{t)  =  Zij  +  t  •  pij.  The 
dependence  of  on  t  is  quadratic,  so  it  is  easy  to  find  a 
parameter  t  =  to  where  Et  has  a  minimum.  We  replace 
^ij  by  Zij  +  to  •  9ij^  but  vertices  which  have  moved  too  far 
(violating  the  constraints)  are  pulled  back.  This  procedure 
is  iterated.  Actually  the  procedure  described  here  is  rather 
unsophisticated.  Improvements  are  the  following:  (i)  If  pij  is 
already  in  a  position  where  it  must  not  move  higher  because  of 
the  boundary  conditions,  but  pij  >  0,  we  let  pij  =  0.  (ii)  the 
same  with  ‘lower’  instead  of  higher,  (iii)  As  is  well  known  in 
multivariate  optimization,  the  direct  gradient  descent  method  is 
usually  not  very  efficient.  We  use  a  conjugate  gradient  method 
for  updating  the  direction  of  descent  in  each  step,  (iv)  In  view 
of  the  large  number  of  variables,  we  use  a  multigrid  method 
for  minimization. 

As  has  already  been  mentioned.  Option  1  leads  to  a  convex 
optimization  problem,  whereas  Option  2  leads  to  a  non-convex 
one.  Numerical  experiments  however  show  that  no  problems 
with  local  minima  occur,  so  the  amount  of  non-convexity 
present  seems  to  be  low.  For  numerical  optimization  in  general, 
see  e.g.,  [25]. 

D.  Measuring  the  deviation 

In  order  to  measure  the  deviation  of  the  smoothed  surface, 
defined  by  vertices  ^ ,  from  the  original  data  pij ,  we  use  the 
tolerance  cylinders  Zij  as  gauge  bodies.  When  scaling  Zij 
with  center  pij,  there  is  a  smallest  factor  Xij  such  that  the 
scaled  cylinder  touches  the  surface  defined  by  the  smoothed 
data  (which  means  that  the  scaled  cylinder  touches  either  the 
x-parallel  or  the  ^-parallel  polylines  defined  by  the  smoothed 
data).  This  factor  Xij  equals  zero,  if  the  vertex  pij  has  not 
moved  at  all,  i.e.,  if  =  pij,  and  the  surface  passes  through 
Zij,  if  Xij  <  1. 

E.  Void  filling 

Missing  data  may  be  treated  in  different  ways.  There  might 
be  a  procedure  for  filling  voids  which  is  appropriate  for  a 
certain  application  and  which  is  applied  before  fair  polylines 
are  used  for  smoothing.  It  is  however  worth  noting  that  fair 
polylines  are  also  capable  of  filling  voids.  If  Zij  is  unknown, 
we  first  attempt  a  crude  guess  at  it  (any  method  for  filling 
voids  is  sufficient  for  that),  and  endow  the  vertex  pij  with 
a  large  tolerance  cylinder  Zij,  so  that  the  vertex  can  move 
(more)  freely  during  optimization. 

V.  Experimental  results 

All  data  sets  described  below  are  conforming  to  the  DTED- 
2  standard,  with  a  post  spacing  of  1  arc  second.  The  horizontal 
and  vertical  error  margins  are  13  and  5  meters,  resp.,  for  all 
vertices.  This  means  that  after  smoothing  the  surface  must  pass 
through  vertically  oriented  tolerance  cylinders  of  base  diameter 
26  m  and  height  10  m  centered  in  the  individual  data  points. 
The  size  of  those  cylinders  here  is  the  same  for  all  points, 
except  for  those  where  no  height  data  are  available.  Such 
points  get  a  large  tolerance  cylinder,  so  that  they  can  move 
freely.  We  implemented  the  smoothing  procedure  in  C++  and 
used  Matlah  to  generate  iso-height  contour  lines. 
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A.  Example  1 

This  data  set  has  a  size  of  800  x  600  posts  and  is  shown  in 
Fig.  7  as  a  shaded  relief.  It  represents  an  area  of  approximately 
24  X  18  km.  The  original  data  set  has  elevation  values  between 
^min  =339  m  and  Zmax  =  1103  m,  and  the  smoothed  data  set 
has  elevation  values  between  Zmin  =  337.5  m  and  Zmax  = 
1098  m.  Details  of  this  example  have  already  been  used  for 
illustration  in  Figures  1,  2,  and  3. 

In  order  to  analyze  the  asymptotic  behavior  of  the  energy 
minimization,  we  run  much  more  iterations  than  necessary  in 
actual  applications.  The  behavior  of  the  network  energy  using 
smoothing  Option  2  is: 


#  iterations 

0 

10 

100 

1000 

Energy 

4817 

727 

628 

623 

Fig.  8  shows  a  histogram  of  vertex  movements  measured  by 
means  of  the  tolerance  cylinders  as  gauge  bodies  (cf.  Section 
IV-D).  The  amount  of  points  whose  movement  is  a  certain 
percentage  of  the  maximum  possible  movement  is  decreasing 
from  left  to  right,  as  only  to  be  expected.  It  is  interesting  to 
note  that  a  substantial  number  of  points  actually  move  as  far 
as  possible  in  order  to  minimize  energy.  This  is  clearly  seen 
by  the  height  of  the  bar  at  the  extreme  right. 

Fig.  9  illustrates  the  movement  of  vertices  during  the 
smoothing  process.  The  height  difference  in  meters  is  coded 
as  a  grey  value.  The  dark  lines  show  valleys,  which  during 
smoothing  are  lifted,  whereas  the  light  areas  show  ridges, 
where  the  smooth  terrain  is  lower  than  the  original  one.  In  both 
cases,  the  deviation  of  the  smoothed  terrain  from  the  original 
one  is  always  consistent  with  the  specified  error/accuracy 
bounds. 

B.  Example  2 

Our  second  example  demonstrates  the  difference  between 
the  two  options  mentioned  in  Section  IV-B.  Option  1  keeps 
the  data  points  inside  the  tolerance  cylinders,  whereas  Option  2 
ensures  that  the  surface  passes  through  the  tolerance  cylinders, 
but  allows  the  vertices  to  leave  the  respective  cylinders 
Zij.  The  results  of  smoothing  are  visualized  in  Fig.  11. 
Minimal  elevation  values  for  the  original  data  set,  and  the  two 
smoothed  data  sets  are  431.00  m,  427.28  m,  and  426.00  m,  re¬ 
spectively.  Maximum  elevation  values  are  724.00  m,  723.27  m, 
and  723.27m. 

While  the  general  effect  of  smoothing  is  clearly  visible 
in  the  contour  plots  for  both  options  (see  Fig.  ll(a)-(c)),  it 
is  hard  to  see  any  difference  between  options  1  and  2.  The 
effect  that  vertices  may  move  outside  their  tolerance  cylinder 
is  strongest  where  the  terrain  is  steep.  In  flatter  areas,  the 
vertices  cannot  move  away  much  from  the  boundary  of  the 
tolerance  cylinder.  This  means  that  the  maximum  deviation  of 
the  smoothed  terrain  from  the  original  occurs  in  those  areas 
where  the  terrain  is  steep.  In  this  example,  this  happens  only 
in  small  parts  of  the  data  set.  Figures  10(a)  and  10(b)  illustrate 
this  deviation,  where  the  absolute  value  of  deviation  is  encoded 
in  grey  values:  Maximum  deviation  is  black,  minimum  is 
white. 

The  fact  that  Option  2  causes  deviation  peaks  in  small  areas 
causes  Fig.  10(b)  to  look  lighter  than  Fig.  10(a),  even  if  the 


Fig.  7.  Data  set  of  Example  1  before  (top)  and  after  (bottom)  smoothing. 


Fig.  8.  Vertex  movement  histogram  for  Example  1. 


Eig.  9.  Vertex  movement  visualization  for  Example  1. 
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Fig.  10.  The  data  set  of  Example  2.  From  left:  Deviation  of  smoothed  data  (Option  I)  from  original;  Deviation  of  smoothed  data  (Option  2)  from  original; 
Difference  between  two  options  of  smoothing;  Vertex  movement  histogram  (Top:  Option  I,  Bottom:  Option  2.). 


absolute  deviation  values  are  greater  in  general.  Fig.  10(c) 
shows  the  difference  between  smoothing  with  Option  1  and 
smoothing  with  Option  2:  the  absolute  value  of  the  height 
difference  is  encoded  in  grey  values  (black  is  maximal  devi¬ 
ation).  The  vertex  movement  statistics  of  Fig.  10(d)  do  not 
show  absolute  height  data,  but  the  deviation  computed  with 
the  tolerance  cylinders  as  gauge  bodies. 

C.  Example  3 

The  third  example  uses  a  small  data  set  of  size  217  x  181, 
which  contains  several  holes  (Fig.  12).  Before  smoothing, 
the  missing  z  coordinates  have  been  guessed  by  a  simple 
averaging  procedure.  During  smoothing,  they  are  allowed  to 
move  freely.  This  example  illustrates  the  fact  that  fair  polyline 
networks  are  capable  of  filling  voids  in  one  pass  together  with 
smoothing  (Fig.  12(b)). 

We  are  now  going  to  further  exploit  the  great  amount  of 
flexibility  which  is  present  in  our  way  of  treating  constraints. 
For  instance,  features  of  the  terrain  which  should  survive  the 
smoothing  process  with  greater  accuracy  than  other  parts  of  the 
data  set  can  be  given  smaller  tolerance  cylinders.  An  example 
of  this  is  given  by  Fig.  12(c),  where  smoothing  has  been 
performed  with  three  different  cylinder  sizes  (diameter/height 
in  meters):  26/10  for  most  points,  2/1  for  the  feature  marked 
in  grey,  and  an  infinite  size  for  voids  in  the  data. 

The  difference  between  the  two  different  ways  of  smoothing 
is  shown  by  Figures  13.  Here  absolute  values  of  z  coordinate 
differences  are  encoded  as  grey  values,  with  white  as  minimum 
and  black  as  maximum.  Data  voids  are  also  shown  in  black. 
The  small  movement  of  vertices  in  the  feature  area  shows  up 
as  white.  A  shaded  relief  of  absolute  value  of  the  difference 
of  the  two  smoothed  data  sets  is  shown  in  Fig.  13(c). 

D.  Example  4 

When  working  with  digital  terrain  elevation  data,  one  might 
want  to  be  able  to  perform  the  following  two  tasks:  (a)  the 
whole  data  set  is  smoothed  with  guaranteed  error  bounds.  On 
demand  a  contour  plot  of  the  smoothed  terrain  is  created,  (b)  a 
single  smooth  contour  line  at  a  certain  fixed  elevation  z  =  H 
is  extracted  from  the  original  terrain  data,  always  respecting 
the  given  error  bounds.  The  first  task  can  be  performed  by 
the  algorithm  described  above  and  is  shown  by  Examples  1- 
3.  The  second  task  can  be  done  by  the  same  algorithm  which 
acts  on  a  subset  of  the  given  data  defined  by  the  inequality 


Fig.  II.  The  data  set  of  Example  2.  From  bottom  to  top:  (a)  Original  data, 
(b)  Smoothing  with  Option  I.  (c)  Smoothing  with  Option  2. 

Hi  <  z  <  H2,  where  [Hi^H2]  is  an  interval  which  contains 
the  height  value  H  under  consideration  and  is  wide  enough  to 
include  the  tolerance  cylinders. 

Figures  14  and  15  illustrate  this  procedure:  We  use  a  data 
set  of  size  436  x  263  containing  both  mountainous  and  fiat 
regions  between  Zmin  =  1894  m  and  Zmax  =  3033  m.  We 
pick  the  2212  m  contour  for  smoothing.  The  region  Hi  < 
z  <  H2  with  Hi  =  2192  m  and  H2  =  2232  m  is  shown  in 
Fig.  14(c).  The  result  of  smoothing  is  shown  in  Fig.  15(a). 
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Fig.  12.  The  data  set  of  Example  3.  From  left:  (a)  Original  data,  (b)  After  filling  voids  and  smoothing,  (c)  After  filling  voids  and  smoothing,  with  features 
preserved. 


Fig.  13.  The  data  set  of  Example  3.  From  left:  (a)  Height  differences  between  original  and  smoothed  data,  (b)  Height  differences  between  original  and 
smoothed  data,  with  features  preserved,  (c)  Difference  between  the  two  smoothed  data  sets  (with  and  without  features  preserved). 


Fig.  15.  The  data  set  of  Example  4.  From  left:  (a)  Smoothing  the  2212  m  contour  by  smoothing  a  neighborhood,  (b)  Smoothed  whole  data  set  with  its 
corresponding  2212  m  contour,  (c)  Deviation  from  original  data  when  smoothing  a  neighborhood  of  a  single  contour  (shaded  relief). 


For  comparison  we  have  also  smoothed  the  whole  data  set 
and  show  its  corresponding  2212  m  contour  in  Fig.  15(b).  In 
Fig.  15(c)  we  show  a  shaded  relief  of  the  deviation  from  the 
original  data  when  smoothing  only  a  neighborhood  of  a  single 
contour. 


VI.  Conclusion 

We  have  presented  an  framework  which  uses  fair  polyline 
networks  for  smoothing  digital  terrain  elevation  data  with 
guaranteed  error  bounds  and  feature  preservation.  The  algo¬ 
rithm  is  capable  of  smoothing  the  terrain  data  with  tolerance 
cylinders  of  different  sizes.  These  flexible  tolerances  have  two 
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advantages  in  particular:  (i)  we  can  preserve  features  present 
in  the  data  by  reducing  the  size  of  the  tolerance  cylinders 
in  feature  areas,  (ii)  the  algorithm  can  be  used  to  fill  holes 
present  in  the  original  data  during  the  smoothing  process. 
Single  contour  lines  are  smoothed  via  processing  of  a  small 
neighborhood  of  that  contour  line. 

Another  important  point  we  would  like  to  make  is  that  the 
smoothing  approach  here  presented  easily  generalizes  from 
height  fields  to  more  general  surfaces,  which  are  becoming 
increasingly  important  in  photogrammetry  and  remote  sensing 
(see  e.g.,  [26]). 
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